
class ewald
{
private:
	double twopioverbox;
	double twopioverbox2;
	double alpha3;
public:
	/*Inicializacion de las sumas de Ewald Coulombiana*/
	void   EwaldCoulomb_init(double length, int particlenumber, double **pos, double **forceR, double **forceK, double *Q, double Ewaldalpha,
			double realspacecutoff, int reciprocalspacecutoff, 	double coulombprefactor, double Qsquare, double epsilon);
	/*Inicializacion de las sumas de Ewald Dipolat*/
	void EwaldDipolar_init(double length, int particlenumber, double **xpos, double **xforceR,
		double **xforceK,double Ewaldalpha, double realspacecutoff, int reciprocalspacecutoff, double epsilon);
	/*Parte directa de las sumas de Ewald de potencial Coulombiano*/
	double EwaldCoulomb_r_space(int force_flag);
	/*Parte reciproca de las sumas de Ewald de potencial Coulombiano*/
	double EwaldCoulomb_k_subset(int start, int end);
	double EwaldCoulomb_k_space(int force_flag);
	/*Componente dipolar de las sumas de Ewald de potencial Coulombiano*/
	double EwaldCoulomb_dipol(int force_flag);
	/*Funcion error complementario*/
	double erfc(double x);
	/*Obtiene el valor de alpha optimo para potencial Coulombiano*/
	void compute_optimal_alpha(void);
	/*Obtiene los qvectors (parte de la interaccion independiente de la posicion de las particulas*/
	void compute_influence_function(void);
	/*Computa las funciones B y C del calculo del potencial Dipolar*/
	void BCFuntions(double r, double *B, double *C);
	/*Computa la funcion D para el calculo de la fuerza del potencial Dipolar*/
	void DFunction(double r, double *D);
	/*Distancia cuadratica y componentes*/
	double distancia2(int i, int j, double *dx, double *dy, double *dz);
	/*Parte directa de las sumas de Ewald para potencial dipolar*/
	double EwaldDipolar_r_space(int force_flag);
	/*Parte reciproca de las sumas de Ewald para potencial dipolar*/
	double EwaldDipolar_k_space(int force_flag);
	/*Funciones de estimacion del valor optimo del parametro de Ewald alpha*/
	double Ewald_cut(void);
	double compute_KPerr_r(double alpha);
	double compute_KPerr_k(double alpha);
	double Ewald_KPerr();
	double Ewald_KPerr_k();
	double Ewald_KPerr_r();
	/*Posicion , fuerza y carga de las particulas (para el caso coulombiano)*/
	double **pos, **forceR, **forceK, *Q;
	/*Numero de particulas*/
	int npart;
	/*Lado de la caja cubica de simulacion*/
	double lbox;
	/*Suma cuadratica de las cargas de las particulas (Coulomb)*/
	double Q2;
	/*Parte real e imaginaria de las exponenciales complejas dependientes de las posiciones de las particulas*/
	double e_re[2*maxkNew+1][2*maxkNew+1][2*maxkNew+1];
	double e_im[2*maxkNew+1][2*maxkNew+1][2*maxkNew+1];
	/*Qvectors*/
	double Ghat[2*maxkNew+1][2*maxkNew+1][2*maxkNew+1];
	double Ewald_get_alpha();
	/*Parametro de corte de Ewald*/
	double  alpha;
	/*cutoff del espacio real y su cuadrado*/
	double  rmax,rmax2;
	/*cutoff en el espacio reciproco y su cuadrado*/
	int     kmax, kmax2;
	/*Prefactor del potencial de Coulomb (Coulomb)*/
	double  prefactor;  /* prefactor of the Coulomb potential */
	double  dipfac;     /* 2*PI / ((1 + 2*epsilon) * L^3) */

	ewald()
	{
	}

	~ewald()
	{
	}

};

void ewald::EwaldDipolar_init(double length, int particlenumber, double **xpos, double **xforceR,
	double **xforceK,double Ewaldalpha, double realspacecutoff,
	int reciprocalspacecutoff, double epsilon)
{
	lbox         = length;
	npart        = particlenumber;
	rmax      = realspacecutoff;
	rmax2     = SQR(rmax);
	kmax      = reciprocalspacecutoff;
	kmax2     = int(SQR(kmax));
	dipfac    = 2.0*PI / ((1.0 + 2.0*epsilon) * lbox*lbox*lbox);
	pos		  = xpos;
	forceR	  = xforceR;
	forceK	  = xforceK;

	alpha     = Ewaldalpha/lbox;
	alpha3=alpha*alpha*alpha;
	twopioverbox=2.0*PI/lbox;
	twopioverbox2=twopioverbox*twopioverbox;

	if (kmax > maxkNew){
		fprintf(stderr,"Error de configuracion de sumas de Ewald:\n"
		"El numero de k-vector requeridi (%d) es mayor\n"
		"que el maximo admitido por el programa (%d)\n"
		"Programa terminado!\n\n",kmax,maxkNew);
		exit(1);
	}

	compute_influence_function();
}

void ewald::EwaldCoulomb_init(double length, int particlenumber, double **xpos, double **xforceR, double **xforceK, double *xQ, double Ewaldalpha, double realspacecutoff,
	int reciprocalspacecutoff,  double coulombprefactor, double Qsquare, double epsilon)
{
	lbox         = length;
	npart        = particlenumber;
	rmax      = realspacecutoff;
	rmax2     = SQR(rmax);
	kmax      = reciprocalspacecutoff;
	kmax2     = int(SQR(kmax));
	prefactor = coulombprefactor;
	Q2        = Qsquare;
	dipfac    = 2.0*PI / ((1.0 + 2.0*epsilon) * lbox*lbox*lbox);
	pos		  = xpos;
	forceR	  = xforceR;
	forceK	  = xforceK;
	Q		  = xQ;

	/* get the alpha value */
	/*if (Ewaldalpha == ALPHA_OPT) compute_optimal_alpha();
	else alpha     = Ewaldalpha/lbox;*/
	alpha     = Ewaldalpha/lbox;
	alpha3=alpha*alpha*alpha;
	twopioverbox=2.0*PI/lbox;
	twopioverbox2=twopioverbox*twopioverbox;

	if (kmax > maxkNew){
		fprintf(stderr,"Error de configuracion de sumas de Ewald:\n"
		"El numero de k-vector requeridi (%d) es mayor\n"
		"que el maximo admitido por el programa (%d)\n"
		"Programa terminado!\n\n",kmax,maxkNew);
		exit(1);
	}

	compute_influence_function();
}

double ewald::compute_KPerr_r(double alpha) {
	/* compute the real space part of the Kolafa Perram error estimate
	*/
	double res;
	res = Q2*sqrt(rmax/(2.0*lbox*lbox*lbox)) * exp(-SQR(alpha)*rmax2) / (SQR(alpha)*rmax2);
	return res;
}

double ewald::compute_KPerr_k(double alpha) {
  /* compute the k space part of the Kolafa Perram error estimate */
  double res;

  res = Q2 * alpha
    * exp(-SQR(PI*kmax/(alpha*lbox)))
    / (SQR(PI) * pow(kmax, 1.5));

  return res;
}

double ewald::Ewald_KPerr() {
  double KPerr_r, KPerr_k;

  /* compute the Kolafa-Perram error */
  KPerr_r = compute_KPerr_r(alpha);
  KPerr_k = compute_KPerr_k(alpha);
  return KPerr_r + KPerr_k;
}

double ewald::Ewald_cut(void)
{
	static double C = -2.837297479480619610825442578061;
	int nx, ny, nz;
	double kx, ky, kz, k2;
	double ghatsum;
	double zeta, zetacutk, zetacutr;
	double Uimgcut;
	double fak1, fak2;

	zeta = C/lbox;

	// compute the sum of the influence function
	ghatsum = 0.0;
	fak1 = 2.0*PI/lbox;
	fak2 = -0.25 / SQR(alpha);
	for (nx = -kmax+1; nx <= kmax; nx++) {
		kx = fak1*nx;
		for (ny = -kmax+1; ny <= kmax; ny++) {
			ky = fak1*ny;
			for (nz = -kmax+1; nz <= kmax; nz++) {
				if (nx != 0 || ny != 0 || nz != 0) {
					if (SQR(nx) + SQR(ny) + SQR(nz) < kmax2) {
						kz = fak1*nz;
						k2 = SQR(kx) + SQR(ky) + SQR(kz);
						ghatsum += 4.0*PI/k2 * exp(fak2 * k2);
					}
				}
			}
		}
	}

	zetacutk = ghatsum/(lbox*lbox*lbox) - 2.0*alpha/wupi;
	zetacutr = 2.0*PI/(lbox*lbox*lbox)*
		(rmax2
		- rmax/(alpha*wupi) * exp(-SQR(alpha)*rmax2)
		- erfc(alpha*rmax) * (rmax2 + 1.0/(2.0*SQR(alpha))));

	Uimgcut = 0.5*Q2 * (zeta - zetacutk - zetacutr);

	return Uimgcut;
}

double ewald::EwaldCoulomb_dipol(int force_flag)
{
	double dip[3];
	double d1=prefactor*dipfac, d2=0.0, aux;
	int    i,k;

	for(i =0;i<3;i++){
		dip[i]=0.0;
	}
	aux=0.0;

	for(k=0;k<3;k++){
		for (i=0; i<npart; i++) {
			dip[k]+=Q[i]*pos[i][k];
		}
		aux+=SQR(dip[k]);
	}

	d2 += d1 * (aux);

	if (force_flag) {
		d1   *= 2.0;
		for(k=0;k<3;k++){
			dip[k] *= d1;
		}
		for (i=0; i<npart; i++) {
			for(k=0;k<3;k++){
				forceK[i][k] -= Q[i] * dip[k];
			}
		}
	}
	return d2;
}

double ewald::EwaldCoulomb_k_space(int force_flag)
{
	int    i;
	double h,c,s,re,im;
	double d1,d2,d3=0.0;


	d1 = prefactor * lbox / (4.0 * PI);

	for (int nx=-kmax; nx<=kmax; nx++){
		for (int ny=-kmax; ny<=kmax; ny++){
			for (int nz=-kmax; nz<=kmax; nz++){
				e_re[nx+kmax][ny+kmax][nz+kmax] = e_im[nx+kmax][ny+kmax][nz+kmax] = 0.0;
			}
		}
	}

	for (int nx=-kmax; nx<=kmax; nx++){
		for (int ny=-kmax; ny<=kmax; ny++){
			for (int nz=-kmax; nz<=kmax; nz++) {
				for (i=0; i<npart; i++){
					if (SQR(nx) + SQR(ny) + SQR(nz) < kmax2) {
						h = 2.0 * PI * (nx*pos[i][0] + ny*pos[i][1] + nz*pos[i][2]) / lbox;
						e_re[nx+kmax][ny+kmax][nz+kmax] += Q[i]*cos(h);
						e_im[nx+kmax][ny+kmax][nz+kmax] -= Q[i]*sin(h);
						//if(cos(h)!=0.0) printf("cos: %20.10le\n",cos(h));
						//if(sin(h)!=0.0) printf("sin: %20.10le\n",sin(h));
					}
				}
					for (i=0; i<npart; i++)	{
						if (SQR(nx) + SQR(ny) + SQR(nz) < kmax2) {
							h = 2.0 * PI * (nx*pos[i][0] + ny*pos[i][1] + nz*pos[i][2]) / lbox;

							c = cos(h);
							s = sin(h);

							d2 = Ghat[abs(nx)][abs(ny)][abs(nz)];
							re = (e_re[nx+kmax][ny+kmax][nz+kmax]*c
								- e_im[nx+kmax][ny+kmax][nz+kmax]*s) * d2;

							if (force_flag) {
								if(nx==0&&ny==0&&nz==0)
								im = (e_re[nx+kmax][ny+kmax][nz+kmax]*s
									+ e_im[nx+kmax][ny+kmax][nz+kmax]*c) * d2;

								d2 = im * Q[i];

								forceK[i][0] += nx * d2;
								forceK[i][1] += ny * d2;
								forceK[i][2] += nz * d2;
							}
							d3 += Q[i] * re;
						}
					}
			}

		}
	}
	d3 *= d1;
			/* self energy correction: */
			d3 -= prefactor * Q2 * alpha / wupi;
	return d3;
}

double ewald::EwaldCoulomb_k_subset(int start, int end)
{
	/* calculates the contribution to the Ewald k-space energy resulting from
	interactions of charges with numbers in {start,start+1,...,end}.
	Note that 'start' and 'end' are INCLUDED in the list! */
	int    i;
	int    nx,ny,nz;
	int    absnx, absny, absnz;
	int    nxpkmax,nypkmax,nzpkmax;
	double k,d1,d2,d3=0.0;
	double Q2_subset=0.0;


	d1 = prefactor * lbox / (4.0 * PI);
	d2 = 2.0 * PI / lbox;

	for (i = start; i <= end; Q2_subset += SQR(Q[i++]));

	for (nx = -kmax; nx <= kmax; nx++)
		for (ny = -kmax; ny <= kmax; ny++)
			for (nz = -kmax; nz <= kmax; nz++)
				e_re[nx+kmax][ny+kmax][nz+kmax] = e_im[nx+kmax][ny+kmax][nz+kmax] = 0.0;

	for (nx = -kmax; nx <= kmax; nx++) {
		absnx = abs(nx);
		nxpkmax = nx + kmax;
		for (ny = -kmax; ny <= kmax; ny++) {
			absny = abs(ny);
			nypkmax = ny + kmax;
			for (nz = -kmax; nz <= kmax; nz++) {
				absnz = abs(nz);
				nzpkmax = nz + kmax;
				for (i = start; i <= end; i++)
					if (SQR(nx) + SQR(ny) + SQR(nx) < kmax2) {
						k = d2 * (nx*pos[i][0] + ny*pos[i][1] + nz*pos[i][2]);
						e_re[nxpkmax][nypkmax][nzpkmax] += Q[i] * cos(k);
						e_im[nxpkmax][nypkmax][nzpkmax] -= Q[i] * sin(k);
					}
					for (i = start; i <= end; i++)
						if (SQR(nx) + SQR(ny) + SQR(nx) < kmax2) {
							k = d2 * (nx*pos[i][0] + ny*pos[i][1] + nz*pos[i][2]);
							d3 += Q[i] * Ghat[absnx][absny][absnz] *
								( e_re[nxpkmax][nypkmax][nzpkmax] * cos(k)
								- e_im[nxpkmax][nypkmax][nzpkmax] * sin(k) );
						}
			}
		}
	}
	d3 *= d1;

	/* self energy correction: */
	double corr = prefactor * Q2_subset * alpha / wupi;
	d3 -= prefactor * Q2_subset * alpha / wupi;

	return d3;
}

double ewald::EwaldCoulomb_r_space(int force_flag)
{
	int    t1,t2;
	double dx,dy,dz,r,r2,ar,ar2;
	double d1,d2,d3=0.0,d4;

	for (t1 = 0; t1 < npart-1; t1++){   /* Quick and Dirty N^2 loop */
		for (t2 = t1 + 1; t2<npart; t2++) {
			r2=distancia2(t1,t2,&dx,&dy,&dz);
			if (r2 <= rmax2) {
				ar2 = SQR(ar = alpha * (r = sqrt(r2)));
				d3 += ( d2 = (d1 = Q[t1] * Q[t2] * prefactor) * erfc(ar) / r );

				if (force_flag) {
					d4  = ( d2 + d1*2.0*alpha*exp(-ar2)/wupi) / r2;
					forceR[t1][0] += d4*dx;
					forceR[t1][1] += d4*dy;
					forceR[t1][2] += d4*dz;

					forceR[t2][0] -= d4*dx;
					forceR[t2][1] -= d4*dy;
					forceR[t2][2] -= d4*dz;
				}
			}
		}
	}
		return d3;
}

double ewald::distancia2(int i, int j, double *dx, double*dy, double *dz)
{
	double ddx,ddy,ddz;
	double r2=0.;

	ddx=pos[i][0]-pos[j][0]; ddx-=dround(ddx/lbox)*lbox;
	ddy=pos[i][1]-pos[j][1]; ddx-=dround(ddy/lbox)*lbox;
	ddz=pos[i][2]-pos[j][2]; ddx-=dround(ddz/lbox)*lbox;
	r2=SQR(ddx)+SQR(ddy)+SQR(ddz);
	*dx=ddx;
	*dy=ddy;
	*dz=ddz;
	if(r2==0)
	{
	    cout<<i<<"  "<<j<<endl;
	    cout<<pos[i][0]<<"  "<<pos[i][1]<<"  "<<pos[i][2]<<endl;
	    cout<<pos[j][0]<<"  "<<pos[j][1]<<"  "<<pos[j][2]<<endl;
	}
	return r2;
}

double ewald::erfc(double x)
{
	double A1 = 0.254829592;
	double A2 = -0.284496736;
	double A3 = 1.421413741;
	double A4 = -1.453152027;
	double A5 = 1.061405429;
	double P  =  0.3275911;
	double T, XSQ, TP;

	T  = 1.0/(1.0+P*x);
	XSQ = x*x;
	TP = T*(A1+T*(A2+T*(A3+T*(A4+T*A5))));
	return TP*exp(-XSQ);
}

 void ewald::compute_optimal_alpha(void) {
	/* use bisectional method to get optimal alpha value */
	double alpha_low, f_low;
	double alpha_high, f_high;
	double alpha_guess, f_guess;

	alpha_low = 0.01;
	alpha_high = 10.0;

	f_low = compute_KPerr_r(alpha_low) - compute_KPerr_k(alpha_low);
	f_high = compute_KPerr_r(alpha_high) - compute_KPerr_k(alpha_high);

	if (f_low*f_high > 0.0) {
		printf("Error: Could not init method to find optimal alpha!\n");
		exit(1);
	}

	do {
		alpha_guess = 0.5 *(alpha_low + alpha_high);
		f_guess = compute_KPerr_r(alpha_guess) - compute_KPerr_k(alpha_guess);
		if (f_low*f_guess < 0.0) alpha_high = alpha_guess;
		else alpha_low = alpha_guess;
	} while (fabs(alpha_low-alpha_high) > ALPHA_OPT_PREC);

	alpha = 0.5 *(alpha_low + alpha_high);
}

 void ewald::compute_influence_function(void)
{
	/* calculates the influence function 2.0/L^2 * exp(-(PI*n/(alpha*L))^2)/n^2
	as a function of lattice vector n (NOT k=2*PI*n/L). This is stored in the
	array Ghat[Maxkmax+1][Maxkmax+1][Maxkmax+1]. For symmetrie reasons only
	one octant is actually needed. */

	int    nx,ny,nz;
	double qua,fak1,fak2;


	fak1 = 2.0/SQR(lbox);
	fak2 = SQR(PI/(alpha*lbox));

	for (nx=0; nx<=kmax; nx++){
		for (ny=0; ny<=kmax; ny++){
			for (nz=0; nz<=kmax; nz++) {
				if ((nx==0) && (ny==0) && (nz==0)) Ghat[nx][ny][nz]=0.0;
				else {
					qua = SQR(nx) + SQR(ny) + SQR(nz);
					Ghat[nx][ny][nz] = fak1 * exp(-fak2*qua) / qua;
				}
			}
		}
	}
}

 void ewald::BCFuntions(double r, double *B, double *C)
{
	double ar = alpha*r;
	double ar2= ar*ar;
	double arpi=2.0*ar/sqrt(PI);
	double r3=r*r*r;
	double r5=r3*r*r;
	double auxB, auxC;

	double erc=erfc(ar);
	auxB=(erc+arpi*exp(-ar2))/r3;
	auxC=(3.0*erc+arpi*(3.+2.*ar2)*exp(-ar2))/r5;
	*B=auxB;
	*C=auxC;
}

 void ewald::DFunction(double r, double *D)
 {
	 double ar = alpha*r;
	 double ar2= ar*ar;
	 double f;
	 f=(15.0*erfc(ar)+2.0*ar/wupi*(15.0+10.0*ar2+4.0*ar2*ar2)*exp(-ar2))/pow(r,7);
	 *D=f;
 }

 double ewald::EwaldDipolar_r_space(int force_flag)
{
	int    p1,p2;
	double dx,dy,dz,r,r2;
	double d3=0.0;
	double B,C,D;
	double coefxy, coefz, Ddz2;
	double fx,fy,fz;

	if (force_flag) {
		/*Inicializamos las fuerzas en espacion directo*/
		for(p1=0;p1<npart;p1++){
			for(p2=0;p2<3;p2++){
				forceR[p1][p2]=0.0;
			}
		}
	}

	/* Quick and Dirty N^2 loop */
	for (p1 = 0; p1 < npart-1; p1++){
		for (p2 = p1 + 1; p2<npart; p2++) {
			r2=distancia2(p1,p2,&dx,&dy,&dz);
			r=sqrt(r2);

			BCFuntions(r, &B, &C);

			if (r2 <= rmax2) {
				d3 += (B-C*dz*dz);
				if (force_flag) {
					DFunction(r, &D);
					Ddz2=D*dz*dz;
					coefxy=C-Ddz2;
					coefz=3.0*C-Ddz2;

					fx=coefxy*dx;
					fy=coefxy*dy;
					fz=coefz*dz;

					forceR[p1][0] += fx;
					forceR[p1][1] += fy;
					forceR[p1][2] += fz;

					forceR[p2][0] -= fx;
					forceR[p2][1] -= fy;
					forceR[p2][2] -= fz;
				}
			}
		}
	}
		return d3;
}

 double ewald::EwaldDipolar_k_space(int force_flag)
{
	int    i,j;
	double h,c,s,re,im;
	double d1,d2,d3=0.0;

	d1 = lbox / (4.0 * PI);


	if (force_flag) {
		/*Inicializamos las fuerzas en espacion reciproco*/
		for(i=0;i<npart;i++){
			for(j=0;j<3;j++){
				forceK[i][j]=0.0;
			}
		}
	}

	for (int nx=-kmax; nx<=kmax; nx++){
		for (int ny=-kmax; ny<=kmax; ny++){
			for (int nz=-kmax; nz<=kmax; nz++){
				e_re[nx+kmax][ny+kmax][nz+kmax] = e_im[nx+kmax][ny+kmax][nz+kmax] = 0.0;
			}
		}
	}

	for (int nx=-kmax; nx<=kmax; nx++){
		for (int ny=-kmax; ny<=kmax; ny++){
			for (int nz=-kmax; nz<=kmax; nz++) {
				for (i=0; i<npart; i++){
					if (SQR(nx) + SQR(ny) + SQR(nz) < kmax2) {
						h = twopioverbox * (nx*pos[i][0] + ny*pos[i][1] + nz*pos[i][2]);
//						e_re[nx+kmax][ny+kmax][nz+kmax] += nz*cos(h);
//						e_im[nx+kmax][ny+kmax][nz+kmax] -= nz*sin(h);
                        e_re[nx+kmax][ny+kmax][nz+kmax] -= nz*sin(h);
						e_im[nx+kmax][ny+kmax][nz+kmax] += nz*cos(h);
					}

	//				 e_re[nx+kmax][ny+kmax][nz+kmax] -= twopioverbox*nz*sin(h);
    //                 e_im[nx+kmax][ny+kmax][nz+kmax] += twopioverbox*nz*cos(h);
				}
				for (i=0; i<npart; i++)	{
					if (SQR(nx) + SQR(ny) + SQR(nz) < kmax2) {
						h =twopioverbox * (nx*pos[i][0] + ny*pos[i][1] + nz*pos[i][2]);

						c = cos(h);
						s = sin(h);

						d2 = Ghat[abs(nx)][abs(ny)][abs(nz)];
//						re = (e_re[nx+kmax][ny+kmax][nz+kmax]*c
//							- e_im[nx+kmax][ny+kmax][nz+kmax]*s) * d2;

                        re = (-e_re[nx+kmax][ny+kmax][nz+kmax]*s
							+ e_im[nx+kmax][ny+kmax][nz+kmax]*c) * d2;

						if (force_flag) {
							im = (e_re[nx+kmax][ny+kmax][nz+kmax]*s
								+ e_im[nx+kmax][ny+kmax][nz+kmax]*c) * d2;

							d2 = im * nz*nz*twopioverbox2;

							forceK[i][0] += nx * d2;
							forceK[i][1] += ny * d2;
							forceK[i][2] += nz * d2;
						}
						d3 +=  nz*re;
					}
				}
			}
		}
	}
	d3 *= d1*twopioverbox2;
	/* self energy correction: */
	double self = 2.*alpha3*npart/3./sqrt(PI);
	//cout<<"value//correccion: "<<d3<<"  "<<self<<endl;
	d3 -=self;
	//cout<<d3<<endl;
	return d3;
}

